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SUMMARY 


A Lagrangian method is developed to solve the Euler equations of gas 
dynamics. The solution of the equations is obtained by a numerical computa- 
tion with the well-known Flux-Corrected-Transport (FCT) numerical method. 
This procedure is modified so that the boundary treatment is accurate and 
relatively simple. Shock waves and other flow discontinuities are captured 
monotonically without any type of fitting procedures. The Lagrangian method 
is employed so that the problem of mesh generation is completely avoided. 

The method is applicable to all Mach numbers except the low subsonic range 
where compressibility effects are small. The method is applied to a one- 
dimensional Riemann problem (shock tube) and to a two-dimensional supersonic 
channel flow with reflecting shock waves. 


I. INTRODUCTION 


One of the outstanding problems in numerical, transonic-flow calcula- 
tions is the exact treatment of boundary conditions. Dramatic changes in the 
transonic flow field are caused by small changes in the boundaries, thus 
exact computation of the boundary conditions is mandatory. For example, a 1% 
change in the thickness ratio of an airfoil can change the total lift coeffi- 
cient more than 10%. Such a strong influence is also seen in the effect a 
body has on the flow over the airfoil in a three-dimensional wing-body con- 
figuration. Furthermore, the curvature at the leading edge of a transonic 
airfoil has a strong influence on the downstream flow field. Thus, accurate 
drag polars, for example, are difficult to compute efficiently by present 
numerical methods. 

The current workhorse for numerical computation of inviscid transonic 
flows is the small-disturbance theory (ref. 1). It can handle steady, 
unsteady, inviscid, and three-dimensional transonic flows. The boundary 
conditions are easily satisfied by the small-disturbance theory, but not all 
types of aerodynamic configurations or flow situations are properly computed 
within the small-disturbance assumptions. The strong feature of the theory 
is the simplified treatment of the boundary conditions. The boundary condi- 
tions (airfoil slopes) are applied on some mean-surface slit instead of the 
actual airfoil surface. Multielement and wing-fuselage configurations are 
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thus easily handled. The main drawback of the small disturbance theory is 
that it is not valid at the wing leading edge, especially wings with blunt 
leading edges, such as the shock-free peaky airfoils (ref. 2), and for wings 
at high angles of attack. Another disadvantage is that the small-disturbance 
shock jump conditions no longer satisfy the Rankine-Hugoniot relations for 
Mach numbers greater than 1.3. 

A more exact approach is to use the full potential equation of transonic 
flow. These include the method of Jameson (ref. 3) and the hodograph method 
of Boerstoel (ref. 2). While both methods can treat the boundary condition 
accurately, they are not straightforward and simple. The Jameson two- 
dimensional procedure requires that the exterior flow field be mapped numer- 
ically into the interior of a circle. The hodograph method is limited to 
two-dimensional flows, with presently little hope for extension to three 
dimensions. Of course the potential equation has the same limitation as the 
small-disturbance equation in that the flows are limited to isentropic flows. 
Thus no shock Mach numbers greater than 1.3 are allowed. 

The Mach number limitation can be overcome by using the full Euler 
equations. This approach is taken by the finite volume method (ref. 4), the 
body coordinate transformation (ref. 5), and the Euler-Lagrange method (ref. 
6). The most promising seem to be the methods of Rizzi and Thompson. How- 
ever, the body coordinate method of Thompson is rather complicated and 
numerically inefficient, especially in the case of an oscillating flap- 
airfoil configuration. The Euler-Lagrange methods of the Los Alamos group 
(ref. 6) are rather disappointing for transonic flows in that the shock 
waves are not captured very well and their schemes require large computa- 
tional effort. The good feature about all the above approaches is that the 
boundary conditions can be included to the same order of accuracy as the 
numerical schemes employed. 

Rather than try to improve any of the above methods, we have chosen a 
different approach altogether. It is a purely Lagrangian approach, unlike 
the methods of Harlow and Amsden, which still use the physical coordinates 
as the independent variables. In the purely Lagrangian approach the inde- 
pendent variables are the Lagrangian coordinates (ref. 7). To obtain the 
Lagrangian equations of gasdynaraics for inviscid, adiabatic and ideal gases, 
the two-dimensional (say) Euler equations of gasdynamics are transformed 
into the Lagrangian divergent free form. To the resulting set of equations 
expressing conservation of mass, momentum, and energy, are added the two 
kinematic relations defining the location (x,y) of a fluid element identified 
by the Lagrangian coordinates (a,b). The numerical computation is carried 
out in the Lagrangian plane. 

The salient feature of this approach is that the physical plane is 
automatically transformed into a rectangular computation plane. No explicit 
transformation computation as in Thompson's method is required. Another fea- 
ture is that the initial distribution of fluid elements or particles is com- 
pletely arbitrary. Thus more particles can be clustered close to the 
aerodynamic configuration for better resolution. The secular mesh distortion 
which usually destroys the accuracy of the Lagrangian calculation does not 
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occur for transonic flows without shear except at stagnation points. All 
these points will be discussed in section II along with the derivation of the 
governing equations in terms of an arbitrarily moving mesh system given in 
Appendix A. 

The main disadvantage of Lagrangian methods are the lack of good numeri- 
cal solution techniques applicable to them. Most of the present "standard" 
numerical methods were developed for Eulerian systems and do not necessarily 
work well for Lagrangian systems. Much experimentation with different types 
of numerical methods was required to arrive at the best method. The results 
of this experimentation are given in section III. 

Several sample problems solved by the present approach are presented in 
section IV. These include a one-dimensional Riemann problem (bursting dia- 
phragm) and the ordinary oblique shock wave reflection off a channel wall for 
a two-dimensional case. 

In the following sections, we restrict ourselves to one or two dimen- 
sions. However, extending the Lagrangian approach to three dimensions is 
simple and adds no extra complexity other than an extra independent variable. 
The three-dimensional Lagrangian equations are given in Appendix B. 


II. FORMULATION 


Governing Equations 

The equations of gasdynamics in two dimensions for an arbitrarily moving 
coordinate system are in conservative form (see appendix A for the derivation), 

U + F + G,_ + H = 0 (1) 
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where 

u,v Cartesian components (x,y) of the fluid velocity 

^m’^m Cartesian components (x,y) of the moving mesh 

a,b coordinates of the moving coordinate system 

x,y coordinates of the fixed Cartesian coordinate system 

^ ~ (^a Yb “ ^b Jacobian of the transformation from the moving 

coordinate system (a,b) to the fixed coordinate system (x,y) 

p, p density, pressure 

p u^ + v^ 

E = rr— + n (polytropic ideal gas) 

(Y - 1)P ^ 

The four equations for x^» x^, ya» and y^ are obtained by differentiating 
the last two equations. They are added to keep the system strongly conserva- 
tive. The first eight equations are sufficient to obtain the velocity, 
density, and pressure field; the last two are required only to locate the 
position of the moving coordinates. 


The fluid velocity components resolved in the moving coordinate system 
(a,b) are given by 

M yb Xb 

ya Xg 

The Jacobian matrices of system (1) are given by 8F/9U and 9G/9U , and the 
system can be rewritten as 



U + AU + BU^ + H = 0 
tab 


( 3 ) 
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with 


A = ^ ; B = ^ 


3U 8U 

The eigenvalues of A are 

, _/M M\ /M M\ / M M ^ c ,/-TT 2^ 

M.2,3.4 - (,'3 - uj. (u - - U^± jVV+ j 


^5, .... 10 - 0 

and those of B are 


t, / M M\ / M M\ / M M ^ c ,/"TT 2\ 

bl, 2 , 3 , 4 =^v -vj, (v - v^), (v - y^) 

^5, . . . , 10 = 0 

where = yp/p and y = ratio of the specific heats. 


The system of equation (1) is the most general set of equations describ- 
ing two-dimensional, inviscid, ideal gasdynamics in any arbitrarily moving 
coordinate system (a,b). The equation describing the conservation of mass, 
momentum, and energy may be written as 


with 



(A) 


-pj - 


" 0 “ 



9 / \ 9 / \ 

pJu 

-► 
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^ VhV - ^ Va^) 

/ V r\ J V 
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-^(V)+ 3 b kO 
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where superscript M is defined by equation (2). The conserved quantities 
are U; (u^ - u^) and (v^ - _v are the (a,b) components of the relative 
convection velocity; and S are the source terms. The eigenvalues of the 
Jacobian matrices of equations (4) are 



6 = 
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It is interesting to note that the eigenvalues depend strictly on the 
convective speed and the relative velocities between the convective speed and 
the sonic speed times some cell distortion factor. This fact remains true 
for any arbitrarily moving coordinate system. Initially the two coordinate 
systems can coincide so that x = a and y = b, but this is not required. 

The equations may now be specialized by specifying the mesh velocity. 

Eulerian system .- The Eulerian system is characterized by a stationary 
M M 

mesh. Setting u = 0 and v =0 obtains x = a, y = b for all time, thus 
m m 

X = 1, y = 1, x^ = 0, y^ = 0, and J = 1. Substituting into equations (1) 

^ o MM 

and (2) obtains (u = u, v = v) 
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The last six equations are trivial and contain no unknowns. The eigenvalues 
of the Jacobian matrices are 


u. 

u 

+ 

C, 

u 

V, 

V 

+ 

C, 

V 


Thus we have recovered the familiar Eulerian description of the two-dimensional 
gasdynamic equations. If we set the mesh velocities to a nonzero constant, 
we obtain a system of equations very similar to the Eulerian description 
except for the convective velocities. Let um = constant, Vj, = constant. 

This results in x^ = 1, y^ = 1, x^ = 0, y^ = 0, and J = 1, since we have 
X = Uj„(t - to) + a and y = Vjj(t - tg) + b. The system of equations becomes 


U 
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Again the last six equations are trivial. The eigenvalues now become 


A = (u - u ) , (u - u ) + c , (u - u ) 
m m — m 

g = (v - V ), (v - V ) + c, (v - V ) 
m m — m 
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We have gained the flexibility of controlling the magnitude and sign of the 
eigenvalues with no extra computational effort beyond the regular Eulerian 
description. However difficulties can arise at the boundaries. 


Parametric system.- Suppose we do not require that initially x = a and 


y = b, and let 


-‘ra» 


Vj[, be arbitrary; a parametric system (ref. 7) is then 


obtained. For comparison with the Eulerian system assume Um = v^ = 0. We 
then obtain 



Jp 


JpuM 


Jpu 


Jpuu^ + yjjP 

u = 

Jpv 

; F = 

Jpvu^ - x^p 
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0 
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The last six equations are the same as in equation (6) and are not shown. 
The eigenvalues of the Jacobian matrices are 


. M M . c /~2~r 2 W 

A = u , u ± j^x2 + y2 , u 


and 


„ M M e / 2 1 2 ^ 

B = V , V ± 3^x2 + y2 , V 


This set of equations is the transformed Euler equations in conservative form, 
where the usual notation is a H ^ and b = n with 5 = C(x,y) and n = n(x,y) 
being arbitrary. 


Lagrangian system.- The Lagrangian description is obtained by setting 
Uji, = u and Vm = V, that is, the mesh moves with the fluid. 
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(Continued) 
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( 8 ) 

(Concluded) 


Here only the first equation is trivial and states that the mass of an element 
of volume Aa Ab remains constant. The eigenvalues of the Jacobian matrices 
are 


M.2,3.4 = 0. ± 3 + x2 . 0 

^ 5,. ..,10 = 0 

and 

61,2,3,4 = 0. ± 3 yfvi + , 0 

65, . . . , 10 = 0 


The Lagrangian description of the gasdynamlc equations in three dimensions is 
similar to the two-dimensional case. It is given in appendix B. 

The eigenvalues of the Lagrangian system are not directly dependent on 
the flow velocities but depend strongly on the mesh configuration. The eigen- 
values are identical, whether the gas is poly tropic or isentropic, and whether 
only the first four equatioijs of the system (8) are used to determine the 
Jacobian matrices (4x4) 3F/8? and 3^/9U with Xa, x^, ya, and yb held con- 
stant, or if the full system of eight equations is used to obtain the (8 x 8) 
Jacobian matrices. 


Before we go into a discussion of the properties and modifications of 
the Lagrangian system (eq. (8)), we will discuss the function of eigenvalues 
in finite difference methods. It is well known that a system of one- 
dimensional hyperbolic equations can be decoupled. If 

Uj. + [F(U)]g = 0 

then 

-► -> 

U + AU = 0 
t a 
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where 



is the Jacobian matrix. We can now find a similarity transformation T which 
will diagonalize A, namely 



The Cj^'s are the eigenvalues of the matrix A and are real for a hyperbolic 
system. The system becomes (T'^U^) + T“^AT(T“^Ua) 0* Assume that T can 
be brought under the differentiation (true, in general, only for the linear- 
ized case) and set u = T~^U to obtain 

u. + Cu =0 
t a 

The system is now decoupled 


(Ui)^ + C^(Ui)^ = 0 

We now solve this system by a Lax-Wendroff scheme. The truncated form 
of the actual differential equation, termed the modified equation by Warming 
and Hyett (ref. 8), being solved by the finite difference scheme is for con- 
stant C 

C 

(ui)j. + Ci(ui)^ + -^(Aa^ - cj At2)(ui)^^^ 

+ ^ - Cl 4t2) + 0(c;(aa2 - C| At2) ] - 0 

where Aa, At are the spatial and temporal step used to discretize the sys- 
tem. The Courant-Friedrich-Lewy number CFL = CAt/Aa < 1 for stability of 
the scheme. If CFL = 1, the numerical scheme solves the original differen- 
tial equation exactly. 

Unfortunately, for a system, the Cj^'s are not constant, thus the 
CFL = 1 condition can be satisfied for only one equation; the others will 
introduce truncation errors. For the Lagrangian system, however, the eigen- 
values can be adjusted by an appropriate mesh configuration. We now choose 
the mesh points so that (the signs of the eigenvalues are Irrelevant) 


""b ^b = f y ''a + ^a - ^ = constant 
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for all points in the computational domain. The modified equation now becomes 


(ui)j. = X(ui)^ + -| (Aa2 - X^At^) (ui)^^^ 


+ (Aa^ - + o|x”(Aa2 - x2At2)j = 0 

Setting the time step to the maximum allowable At = Aa/X results in 

(ui)^ + ^(ui)^ = 0 


In other words, the finite difference scheme solves the system of partial 
differential equations exactly. Two sources of error can still occur; namely, 
the inaccuracies in determining the mesh configuration and nonlinearity. The 
above analysis is for the linear case only. For the nonlinear case such as 
the invlscld gasdynamic equations, the above procedure of mesh adjustment 
will increase the accuracy of the numerical scheme by at least one order. 

This can be seen from Lerat and Peyret’s analysis of the numerical solution 
of the inviscid Burger's equation by a generalized Lax-Wendroff scheme. 

The above ideas are still in flux and no numerical verifications have 
been made. This is the end of the digression and now some properties and 
variations of the Lagrangian system (8) will be given. 

Jump conditions for Lagrangian system .- Let o(a,b,t) = 0 define a dis- 
continuity surface in the Lagrangian space, then the jump conditions of system 
(8) are 

'“1 1? ^ Is 1^1 f ■ ° 


where [f] indicates the jump in the value of f in crossing the discontinuity 
surface a = 0. The last two equations of system (8) are not included since 
they do not enter in the solution of the previous eight equations. 


Total enthalpy form of the energy equation .- Consider the energy equa- 
tion 


or setting on 


[PJE] If + [JpuMj If + [JpvM] If . 0 
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u 



+ V 
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j^(pJE 1^ + JpqM • Va)]= 0 
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We have 


and 


9g _ _9a ^ 
9t ~ 9a 9t 


+ i£ Ik 

9b 9t 


da 

dt 


9a . 9a . 9a 
o 9x a 9y 9t 


hence 


li 

9t 


9a 9a 

= - u -r V — 

o 9x o 9y 


and similarly for 


9t 


9b ^ 

9x a 9y 


The Cartesian velocity components of the discontinuity surface in the 
Lagrangian coordinate system are Uq and Vg. Substituting into 9a/9t re- 
sults in 


£ / yb% ~ yg \ 90 / -ya^g ^a^ \ 9g 
t ~ \ J /9a“\ J /9b 


k£ 

9t 


or 


9g M 9£ M 9£ 
9t ■ 9a ” "^0 9b 


'M 


= -q • Vo 
^ o 


( 10 ) 


M M 

where u^j and Vg are the Lagrangian (moving mesh) components of the discon- 
tinuity surface in the Lagrangian coordinate system. The discontinuity veloc- 
ity in the Lagrangian system can be expressed in terms of Cartesian velocities, 
that is, 

-♦■M -*■ M ->-M 

q = q - q 

^o ^oc ^ 

->■ M ->M 

where qgg> q are the discontinuity and mesh velocity, respectively, in the 
Cartesian system given, however, in terms of the Lagrangian components. The 
mesh velocity is, of course, the fluid velocity. Hence 


9o 

= q • Vg - q 
9t ^ ^oc 


Vo 


The energy jump condition is now 
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( 11 ) 


f)lf = 0 

In terms of the differential equation this becomes 

i Jp(e + ^)+ -^(jp + ^(jp = 0 

If the flow is stationary, =0, the energy equation reduces to 


( 12 ) 


(13) 


The equation shows that the total enthalpy h^ = E + p/p is constant along a 
particle path if the flow is stationary. Since equation (13) was derived 
using the proper jump conditions (eq. (9)), we have the condition that the 
total enthalpy is not discontinuous across o = 0, that is, 

[htl = 0 (14) 

for steady flows, even for Lagrangian systems which always contain time deriv- 
atives . 

Pressure form of energy equation .- The energy equation of system (8) 
can be reduced to 

^0"p) ' “ 


by combining the equation of system (8), assuming the polytropic (not isen- 
tropic) gas relation 


E = 


(Y- Dp 


+ V^) 


and noting that 


at 



^ (-ypu 0 


Another form of equation (15) is obtained by combining it with the continuity 
equation and the state equation relating the entropy in terms of the pressure 
and density for polytropic gases 


to obtain 


s 


+ C, 



^ (pJs) = 0 


(16) 
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The physical meaning of the energy equation (15) is that the entropy, s, of a 
fluid particle is a constant for adiabatic flows. This is as expected if the 
flow is continuous. Unfortunately, equations (15) or (16) do not satisfy the 
jump conditions (eq. (9)). The proper form of equation (16) must be derived 
from the entropy transport equation which is not included in the Euler Equa- 
tions of gasdynamics. The Lagrangian form of the entropy equation for adia- 
batic flows is 

(pJs) > 0 (17) 

which has the jump condition 

[pJs]||>0 (18) 

which is quite different from that of equation (16). However, equation (15) 
or (16) is still useful for transonic flows where the entropy changes are of 

0(Ap3). 


Lagrangian Coordinates 

The Lagrangian coordinates serve to Identify the fluid particles. At 
some initial time a distribution of particles (or cells) is chosen in the 
physical plane (fig. 1). Each particle is identified by the Lagrangian co- 
ordinates (aj^,bj). For simplicity we set a± = i, and bj = j where 
i = 1,2,...,L and j = 1,2,...,M with the total number of cells being LxM. 
The Lagrangian plane will then consist of unit square cells with any bound- 
aries appearing as slits as in figure 1. Single or multielement boundaries 
are possible. Unlike the procedure followed by Lamb (ref. 7), the initial 
distribution of points is not uniform, that is, a kjx, b ^ k 2 y, ki, k 2 
are constants. This nonuniform distribution allows the particles to be 
clustered around boundaries and be identified by integer values, which sim- 
plify the numerics . 


Boundary Conditions 

There are four types of boundary conditions to be considered: inflow, 

outflow, free stream or far field, and impermeable surface condition. All 
boundary conditions are first predicted by the FCT method and then the follow- 
ing procedures are utilized to correct the predicted values. This is very 
similar to Abbett's procedure (ref. 9), except his is restricted to steady 
supersonic flows. To predict the boundary terms by the FCT method without 
turning the flux limiter off, the flow field is linearly extrapolated two 
points beyond the boundary. The procedure will be derived for two dimensions 
but it is also valid for three-dimensional flows. 

The inflow boundary condition is simple. The dependent variables of 
system (8) are set to the specified inflow conditions. The outflow boundary 
is not so simple. If the flow is supersonic, then the particles along with 
their properties are allowed to flow out undisturbed, that is, a linear 
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extrapolation will do and will not disturb the upstream conditions. If the 
outflow is subsonic, then some outflow boundary conditions are required to 
set the dependent variables of system (8). What these outflow boundary con- 
ditions should be is yet unknown. 

The far-fleld boundary condition will be treated similarly to the imper- 
meable surface boundary and this requires that some analytic representation 
of the far-field streamline be given. Let both the far-field streamline and 
the surface streamline be given by 


y = f(x) (19a) 

We have used the term streamline loosely. It includes particle paths and 
thus the streamlines coincide with the Lagrangian coordinate b = constant. 
Differentiation (eq. (19a)) with respect to time obtains 



y^ = f'(x)xt 

(19b) 

We assume that f(x) is differentiable as many times as necessary, 
last two equations of system (8) we have 

From the 


X 

rt 

II 

C 

(19c) 


Yt = V 

(19d) 

thus 

V = f ' (x)u 

(19e) 

Differentiating equation 
respect to "t" obtains 

(19a) with respect to "a” and equation (19b) 

with 


Ya = f '(x) • Xa 

(20a) 


v^ = f ' (x) • Uj. + f"(x) • u2 

(20b) 


The two momentum equations of system (8) are cross multiplied and summed to 
eliminate the p^ terms to obtain 


(pJu)^ 


pJ / 

1 + f'2 ^x^p 


+ u2 f ' 



=- 0 


where use has been made of equations (20a) and (20b). The time may now be 
eliminated by the x-momentum equation of system (8) to obtain 


x^ (1 + f ' 2) - (Xjj + y^^ f ' ) Pg + (pJ) u2 f " = 0 


( 21 ) 
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This is the normal momentvim equation in the Lagrangian coordinates. Since 
we are solving system (8) by a time-splitting procedure, the boundary values 
along the two boundary streamlines are required only during the "b" sweep. 
Hence we are only solving part of system (8); namely. 


1 ^ 1 ^- 


0 


Only G is required along the two streamlines. G is determined by u, v, x , 
ya, and p. We have equations (19a), (19b), (20a), and (21) available. To ^ 
obtain the remaining relation we need to specialize the equations to steady 
or unsteady and isentropic or poly tropic flows. 

Polytropic unsteady .- If the flow is unsteady and polytropic, then we 
simply extrapolate x and u to surface, y is computed and Xg, x^ and y^ are 
obtained by second order accurate finite differencing, either central or 
one-sided as required. The p^ of equation (21) is approximated by a one- 
sided difference and the pg by a central difference to obtain a scalar 
tridiagonal system for the unknown p at the surface. The density is 
obtained from the continuity equation 


p (Pi) 

(Wb - *bya) 

V from equation (19e) and E from the equation of state. 

Polytropic steady .- If the flow is steady, then use is made of the 
steady-state version of the energy equation 


= YP . u^(l + f'^) 
t " (y - l)p 2 


constant 


( 22 ) 


Substituting this into the normal momentum equation to eliminate the u^ 
term obtains 


O 2vJf" 2(pj)f"h^ 

Xg(l + f'2)pb - (Xb + yb f’) pa - , i)(i + f»2) P = (i + f’2) (23) 


The same procedure may now be followed as previously except u is not 
extrapolated. Once p, Xg, yg, Xb, and yb have been obtained, u and v 
result from equations (22) and (19d), respectively. The determination of 
E is also the same as previously described. 

Transonic unsteady or steady .- The only additional simplification in 
transonic flows is that it can be assumed to be isentropic, thus 

p = kp^ (24) 
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The density is now computed from this relation and is computed from 


y 


b 


to keep the system consistent. 




Parametric form of y = f(x) .~ If the bounding streamline behaves such 
that f and/or f' do not exist, then equations (19a, b) and (20a, b) are 
replaced by their parametric representation. Let the arc length along the 
curve y = f(x) be given by a, then equations (19a) and (19b) become 


y = n(a) \ 
X = C(o) I 


vx = u*y 
a a 


and equations (20a) and 20b) become 


y X = X y 
a a a a 


(19a') 

(19b’) 


(20a’) 


x^ V 

a t 


= y x'^ + 

oat 


(x y - y X ) 
a aa a aa 


(20b’) 


Equation (19a’) is the parametric counterpart to equation (19a) and n and 5 
are given continuous and smooth functions of a. We require that the function 
be smooth so that y and x remain finite. The continuity condition guar- 
antees the boundedness of y^j^and x . These restrictions mean that the 
streamline is continuous without any kinks or sharp corners. 


III. NUMERICAL SOLUTION METHODS 


The numerical solution of system (8) is obtained by an explicit hyper- 
bolic differential equation solver. Of the ten equations of system (8) only 
nine need by solved by the D.E. solver. Unfortunately, most of the present 
numerical schemes have been devised for Eulerian systems and it does not 
necessarily follow that they will be efficient for Lagrangian systems. The 
procedures devised for Lagrangian or Euler-Lagrangian systems, such as by the 
Los Alamos group (ref. 6) are not applicable to system (8). Most of the vari- 
ous Lagrangian procedures devised by the Los Alamos group and others operate 
on the Cartesian space and not the Lagrangian space as required by system (8). 
We also required that any discontinuities in the flow field be captured accu- 
rately without excessive and spurious nonphysical oscillations by the numeri- 
cal procedure, as we did not want to add the additional complications of shock 
fitting procedures. 

We finally required that a time-step splitting procedure be used so that 
the number of equations to be solved during each sweep in the two spatial 
directions is reduced to six. This requirement is necessary to make the 
numerical computation in the Lagrangian system competitive with that in the 
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Eulerian system. The Eulerian system requires that four equations be solved 
during each sweep. For these reasons several different numerical schemes 
were tried of which only four will be given here. They are the IlacCormack, 
generalized Lax-Wendrof f , Shasta flux corrected transport FCT, and the "fourth- 
order phase error" FCT schemes. All schemes_^will be given for the one- 
dimensional, nonlinear, hyperbolic system f^- + [F(?)]x+ S = 0. Here 
f(x,^) _^is a column vector of m components and the Jacobian matrix 
A = F'(f) has m real and distinct eigenvalues A(f). 


MacCormack and Generalized Lax-lJendroff Schemes 


The MacCormack schemes are special cases of the generalized Lax-Wendroff 
schemes devised by Lerat and Peyret (ref. 10), thus only the latter schemes 
will be given. The schemes are given in a predictor-corrector form: 


P: 


C: 


where 



(1 - B) f^" + 6 f - ao (F - f" ) - oAtS 


0 

2a 

[(a - B) + 

(23 - 

1) 


X n 



■ a 

- « "i-i ^ - 

^i-1 

] ■ 2a 



( 1 ) 


t = nAt 
X = iAx 
a = At /Ax 


a,g = two arbitrary parameters with a ^ 0 

The predictor f^ is an approximation to the solution at x = (i + 3) Ax and 
t = (n + a) At. The schemes are three-point, second-order accurate for all 
values of a 0 and 3, and stable for ^ 1* Setting a = 1 and 

3 = 0, 1 yields the two MacCormack schemes. Lerat and Peyret determined that 
the optimum values are a = 1 + /5/2(=2.118) and 3=1/2 using the inviscid 
Burgers equation as the model system in their analysis of the truncation 
errors. They were not able to analyze the truncation errors using the full 
Euler equations. If the scheme (1) is applied to a linear system fj- + cf^^ = 0, 
then the actual differential equation, called modified equation by Warming 
and Hyett (ref. 8), solved by the scheme is 

f + cf + ■§• (Ax^ - c^At^) f + ^ (Ax^ - c^At^) f + ... = 0 

t X 6 XXX 8 xxxx 


Nothing can be said about the parameters a, 3 from linear analysis. This is 
not the case, however, for the nonlinear system as shown by Lerat and Peyret. 
For our Lagrangian system, both the MacCormack schemes and the optimum Lerat 
and Peyret scheme were tried with equally disappointing results to be shown 
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in section IV. The introduction of an artificial viscosity to the scheme (1) 
did not improve the results, in fact, the results were even slightly more 
degraded. 


Shasta Flux Corrected Transport Algorithm 


The flux corrected transport algorithms developed by Boris and Book 
(ref. 11) require that the nonlinear hyperbolic system be written in the form 

fj. + (uf)^ + S = 0 

-► 

where u is the convective velocity transporting f. It is similar to the 
previous form of the hyperbolic system. Thus we have for the one-dimensional 
Euler Equation 



p 


0 



— >• 


f = 

pu 

; and s = 

Px 


e 


_(pu)^_ 


where for a polytropic ideal gas 


e 


P pu^ 

Y-1 2 


The Shasta FCT scheme (ref. 12) is as follows. Let f^, u^, and S^ 
be given for i = 1,2, ..., N. First, the untransported and a diffused 
solution are computed 


f 


0 

i+i 


1 






where the superscripts given here and later vzill stand for 

0 - initial 
D - diffused 
T - transported 
TD - diffused and transported 
C - corrected 


Now determine the convective transport factors 
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and source terms for S = 9p/9x 


= o(pl^^ - pj) 


or If S = p, then 


’i+i = f-(Pi+i + pS) 


The transported diffused solution is 


tf - 4(Qt)^ - 4 (qI)^ 


+ Qtwj - Sj+j) + Qi(f;; - s..j) 

The changes due to transport alone are 


T TD D 

6 f! = f;: - f. 

1 i 1 


and differences in the diffused transported fluxes are 

. ,TD -TD TD 
^ i+i ■ ^i+1 ■ ‘i 

and initial fluxes are modified by the convective transport 


These fluxes are now corrected by the nonlinear flux limiter 

= Sgn • Max jo, Min^Sgn • 6 I ^i+i ^ ^ ^i+|)] 


where 

Sgn = sign(f^^^) 


The final flux corrected transported solution at the new time step is obtained 
by 


f 


1 

i 





The corrected flux 
TD 


xcTD -T 

'i+i> *1+1 


Q 

f^^ has four different possible values; namely, 0, 


Therefore, the final time corrected transport solution 

has a total of 16 different finite difference representations. We will ana- 
lyze only the two most prevalent ones; namely, the cases for 
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0 


±+i h-i 


and 




= 

±+i ’ 1-i i-i 


ORIGINAL PAGE IS 
OF POOR QUALITY 


To analyze the Shasta FCT scheme we will drop the source terms and assume 
the Uj^'s to be constant = c. Define a = uo. For the first case (f^+i = 0) 
we have 

f ■ " f™ 
i 1 


- 'i-i) 

which has the modified equation 


e . f qAx 
^t 8 4 


+ c / Ax^ 


1 


c^At^^ f + 

* XXX 




f + . . . = 0 

XXXX 


The modified equation shows that for this case, that is, null flux correction, 
the scheme is only first-order accurate. For the second case 
we have 


- ^U) 

r /f? - f9 A oAfO , - I6f9 T + 

• n - ^ 

T [-1 ^ f?-i) + 1 


16f 


i+1 


- 2f 


i+2 


16f 


i+1 


- 30f° + 16f9 , - f 

i x-1 


12 


«] 


The modified equation is 


ffc + cf^ + 


c^At 

8 


f (¥ - ‘ 

^Ax^ - c^At^^fjj 


XXX 


+ . . . = 0 


Here the scheme has second-order accuracy and this modified equation shows 
that the dispersion error is approximately only one-fourth the dispersion of 


the Lax-Wendroff scheme as shown *by the coefficients of the f 


XXX 


term of 
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the modified equations of the two schemes. This is one of the reasons for 
the improved performance of the FCT scheme over the generalized Lax-Wendroff 
schemes. The standard Fourier (Von Neumann) stability analysis of the two 
methods confirms this conclusion. Figure 2 shows the amplification factor 
I? I and dispersion (or phase) error <j>/4>e the two methods at two CFL 

numbers. For both CFL numbers shown the dispersion error is substantially 
less for the FCT method. Another reason is that the flux limiter sets the 
flux correction fi+^ to zero near discontinuities and thus insures that the 
solution remains monotonic near these discontinuities, where the Lax-Wendroff 
schemes usually show spurious nonphysical oscillations. The FCT scheme has 
only first-order accuracy in this situation. However, this condition usually 
occurs only in a small neighborhood of the discontinuity, so the FCT scheme 
is essentially second-order accurate over most of the computational domain. 

The co^ition for stability is julmaxO 1 1 and for monotonicity is 

I n Imax ° - *^3/2 . 

Although it is not immediately apparent from the scheme given, the essen- 
tial features of the FCT scheme are: 

1. The conserved quantities f^ are changed by the convective trans- 
port and source terms. Enough nonphysical damping is provided so that the 
solution remains monotonlc if initially monotonic. 

T 

2. Flux antidiffusion terms fi+i are computed so that the excessive 
damping of step one can be eliminated where not needed. 

3. However, before the^solution is antidiffused, the antidiffusion terms 
fj^+^ must be modified to fi+J by the nonlinear flux limiter to preserve the 
monotonicity of the antidiffused solution of step one. 

The effect of these three steps of the Shasta FCT scheme is shown graph- 
ically for a one-dimensional Riemann bursting diaphragm problem in figures 3a, 
b, c. These figures show the pressure profile at the initial state and at a 
later time. The first figure shows the solution after the diffused transport 
step. This solution is smooth and monotonic. The antidiffused solution with- 
out the flux limiter is shown in figure 3b. Large overshoots occur near the 
shock which has been steepened from that of the first step. The last figure 
shows the final solution where the antidiffusion fluxes have been limited by 
the nonlinear flux limiter. The solution is nearly monotonic and most of the 
shock jump occurs over two cells compared to six or seven cells after the 
first step. The final solution is not quite monotonic, and that is due to 
the pressure being determined from the three conserved quantities: density, 

momentum, and total energy, all three of which are computed monotonically by 
the FCT scheme. The three phases of the FCT algorithm are more obvious in 
the third order FCT described below. 


Third-Order FCT 

The third-order FCT scheme is the Latest and most accurate scheme devised 
by Boris and Book (ref. 11) who call it the "Fourth-Order Phase Error" FCT. 
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The procedure shoim here differs from that of Boris and Book in two respects: 
(1) the computational mesh is rectangular and fixed but not necessarily uni- 
formly spaced, and (2) the cell node positions are shifted by 1/2 cell width 
from that of Boris and Book. The notation of the mesh configuration is shown 
in figure 4. 


Figure 4 has N grid points and N - 1 interfaces. The interface posi- 
tions are denoted by 



To determine the interface fluxes, we need the interface velocities 



and the interface conserved quantity 


The transported 



i = 1, . . . , N - 1 


fT 

1 


due to convection and source terms can now be computed 


A» = A« f» - 

for i=2, ...,N-1. And for i = 1, N 

fl^= f? 

1 1 

We assume that the source S is given by 8p/8x. Extensions to other types 
of source terms are obvious. The cell volume is 



The diffused transport is now determined by 


TD T 

A.f;|^ = A. ft + v.^A.j, 
1 i XI i+i i+i 




TD 0 

for i = 2, ..., N - 1, and set f^^ = f^ for i = 1, N. The interface vol- 

umes are defined as 


^i+i = - r^ for i = 1 N - 1 

The diffusion coefficient will be defined later. The uncorrected antidiffu- 
sive fluxes are defined as 
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i+i '■i+1 


(^i+1 ■ ^ for i = 1, ..., N - 1 


where u | i is the antidiffusion coefficient. The raw antidif f usive fluxes 
are now corrected. 


^ ^ = Sgn • Max \ 


where 


{o. Ml„ . Sgn . ' f^i). Sgn • A^(ff - f 


(Cx - ^T) 


Sgn = sign 


for i = 2, ...,N-2. The corrected fluxes F 3/2 ^re determined 

by the above relation, dropping the undefined terms outside the boundary. 

The final results are now obtained 


i - ^i Ai [^±+i - ^i-i ) 


for i = 2, ..., N - 1. The diffusive and antidif fusive coefficients are 

2 

- 1 ... 

'^l+i 6 3 


with 


1 €l+i 

^i+i 6 " 6 


. _ r ^t 

^i+i ^ “i+ i 2 




These choices make the scheme third order accurate as shown below. 


If we now assume that we have uniform mesh spacing so that Ax = Ar = 
l/2(ri+3 - r^-x) = and constant convective velocity ui = C, then the 

modified equation of the scheme is 



where e is the CFL number = CAt/Ax. With the above choices of v and y 
the modified equation becomes 
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Ax'* 

+ cf + ^ 4r(«^ - l)f + 
t X 24 At xxxx 


0 


4 

That is, the scheme is third order accurate, 0(Ax /At). In the neighborhoods 
of discontinuities the flux limiter essentially sets p = 0 and the scheme 
becomes accurate of 0(Ax^/At). It can be shown that the CFL restriction for 
stability is CFL < 1 and for monotonicity also <1. 

For nonlinear problems, both of the FCT algorithms require the solution 
to be advanced in two half-time steps to arrive at the final time (n + l)t 
where At < Ax/Amav to maintain the same order accuracy for nonlinear prob- 
lems as for linear ones. This procedure is equivalent to the predictor- 
corrector steps of the Lax-Wendroff schemes. 

Another comment is that for the Lagrangian system (8) there are no con- 
vective velocities and thus it might be argued that the FCT algorithms are 
not applicable for this system. However, the convective transport can be 
considered to be by the eigenvalues of the Jacobian matrix rather than the 
actual fluid convective velocity. For example, Boris and Book. (ref. 11) argue 
that for a one-dimensional Eulerian system, the effective convective veloci- 
ties are u, u, and (1 + p/e)u for the continuity, momentum, and energy 
equations, respectively. Thus, even for a Lagrangian system, the energy 
transport is by the convective velocity (p/e)u. It is more valid, however, 
to consider the eigenvalues to be the effective transport velocities. Con- 
sider the Eulerian system in one dimension 


where 


f 


t 


+ Wx = 0 



w 


pu 

pu^ + p 
(e + p)u 


Now define the Jacobian matrix A = 9w/8f, to obtain 

fj. + (Af)^ = 0 


where we have used the homogeneous property of the Eulerian system (ref. 13). 
A similarity transformation T can now bedefined so that Jacobian matrix A 
is diagonalized and f = T~^f where the fi's differ from the fi's only 
by a constant factor, that is 


D 



0 

u+ c 
0 


0 

0 


u- c 


and f = kif^. The system after some manipulation becomes 
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+ (DT-lf)x - (Tt^f + Tx^w) = 0 
or 

+ (Df)x + S = 0 


where S are the source terms 

S = + T-^w) 


ORIGINAL Page is 
OP POOR QUAUTy 


The inverse similarity transformation is 


1 


y - 1 







uc 

Y - 1 


u 



p 


With this transformation we can compute f to be f = T“^ f H f, so the con- 
stant coefficients are unity. The source terms are 


S = 



_ 3 _ 

3x 




- e [ c + (y 



Here we have a system, completely equivalent to the original system, which 
gives the proper jump conditions and which has the eigenvalues as the effec- 
tive convective velocities. This shows that the convective terms of the 
Eulerian system play no major role in the FCT algorithms. Similarly for the 
Lagrangian systems, which have some nonzero eigenvalues, but zero convective 
velocities. The lack of convective velocities does not prohibit the use of 
the FCT algorithms. 


IV. NUMERICAL RESULTS 


In this section the numerical methods of section II are compared for a 
one-dimensional Riemann problem (bursting diaphragm). The most effective 
method will then be applied to a two-dimensional problem of the ordinary 
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oblique shock wave reflection off a channel wall. To check the boundary 
conditions, the Riemann problem consists of a shock tube with reflecting end 
walls. For the two-dimensional problem, the boundary procedures described in 
section II are used even for the case of straight walls where reflection 
principles would be much simpler. 

The initial conditions and geometry of the Riemann problem are shown in 
figure 5. As a standard of comparison, this problem was solved in the Eulerian 
coordinates with the MacCormack scheme. Initially the pressure profile is a 
step function as in figure 6. After the burst a shock wave propagates to the 
right and a rarefaction wave to the left. At later times these two waves 
reflect off the end walls. The shock and rarefaction waves are resolved quite 
well except for the slight spike of the reflected shock wave and slight oscil- 
lation at the foot of the rarefaction wave. The spatial step Ax = 0.002 is 
quite small and good resolution is to be expected. Changing to the Lagranglan 
system with everything else fixed results in the solution shown in figure 7. 

The oscillations at the root of the expansion wave are larger and the spike 
at the reflected shock wave is much larger. The particles are at first 
equally spaced at Ax = 0.002 so the masses in each cell (pJ = pxg) are 
not uniform initially. The temperature (p/p) profiles are shown in figure 
7b. 


If the initial particle distribution is such that we have equal mass in 
all the cells, then the results shown in figure 8 are obtained. Large oscil- 
lations trail from the shock waves. The expansion waves seem to be resolved 
better. Thus nothing is gained by adjusting the fluid cells to contain equal 
masses . 

Relaxing the restriction of equal masses and using the generalized Lax- 
Wendroff method with the optimum values of the two arbitrary parameters of 
Lerat and Peyret (ref. 10) obtains a much improved solution (fig. 9). The 
oscillations at the shocks are much smaller, but the expansion wave oscilla- 
tions are about the same as previously. Including an artificial dissipation 
of the Lapidus type (ref. 10) does not improve the solution as it does for an 
Eulerian calculation (fig. 10). Other values of the two parameters did not 
improve the results. 

Somewhat improved results are obtained from the FCT method of Book, Boris 
and Hain (ref. 12). Figure 11 shows that the oscillations at the expansion 
wave and outgoing shock wave are gone. A spike appears at the reflected shock 
and also at the contact surface. The contact surface oscillations disappear 
if equal-mass fluid cells are chosen. The spike appears at the reflected 
shock since the flux limiter of the original FCT method is partially turned 
off at the boundaries. 

All the methods tested tend to give comparable results with the FCT 
method giving slightly better results. However, all the cases considered so 
far have a very fine mesh distribution of Ax = 0.002 or 500 points. For 
two or three dimensions this is excessive. Reducing the number of points to 
50 obtains quite different results. These results are shown in figures 12 
through 15 for the optimum generalized Lax-Wendroff method for the Eulerian 
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coordinates and Lagrangian coordinates, and for the 1/8 (Shasta) FCT, and 
third order FCT method with Lagrangian coordinates, respectively. The re- 
sults for the Lax-Wendroff method are not very satisfactory, but for the 
Eulerian calculation additional damping can be added to the scheme to control 
the large oscillations. This does not work for the Lagrangian computation. 

The FCT calculations are shown in figures 14 and 15. Figure 15 gives 
the results of the third order FCT method using Boris and Book's (ref. 11) 
version of the Lagrangian code. This code keeps the cartesian coordinates as 
the independent variables. The FCT results are much better than the Lax- 
Wendroff results. The third order FCT is slightly better than the Shasta FCT 
except for the terracing of the expansion wave. This is because Boris' form- 
ulation does not compute the x^, xj,, ya» Yb explicitly as in the present 
method. 

The second problem solved was a two-dimensional reflection of an incident 
shock wave off a straight wall. The initial conditions and geometry are pre- 
sented in figure 16. The incident shock wave is generated by the 1:5 wedge 
placed along the bottom wall of a straight channel. This shock wave reflects 
off the top wall as shown. The flow properties are chosen so that the reflec- 
tion is regular and so that the flow is supersonic throughout thus the bound- 
ary computation of the present method can be compared with Abbett's scheme 
(ref. 9). Abbett's scheme is valid only for steady supersonic flow, so only 
the steady state results are shown. Another reason for choosing supersonic 
flow is that then the outflow conditions do not influence the upstream flow 
field . 

The initial condition for computation is simply a straight channel with 
a Mach number 2.2 flow throughout. At t = O’*" the bottom wall aft of x = 1 
rotates upward about x = 1 until the final ramp angle of tan“^ (1/5) is 
reached. The computation proceeds until all the transients are washed down- 
stream, which at this Mach number required a particle travel of approximately 
2.2 ramp lengths, after the final ramp angle was attained. In physical time 
the ramp growth required 1 sec and the steady state was reached in another 
second. The total CPU time on a CDC 7600 for convergence was 135 sec, and 
the total number of time steps amounted to 300. 

The results in terms of pressure profiles are shown in figures 17 and 18, 
for the Prandtl-Meyer (Abbett's scheme) and the normal momentum equation 
boundary correction. The oblique shock does not reach the foot of the ramp 
since a fairing function was fitted between x = 0.9 and 1.1 to eliminate 
the discontinuity in the surface slope and curvature at x = 1. The fairing 
function is a quintic polynomial whose value, slope, and curvature matches 
the prescribed boundary at the fairing positions. 

The shocks are captured quite well, usually over three or four mesh 
points. The total number of mesh points is 51 x 21. It should be mentioned 
that the shocks cross the cells at an angle of approximately 45°, so there is 
no mesh alignment with the shock wave. The Prandtl-Meyer boundary correction 
gives slightly better results than the normal momentum correction; the latter 
has slight overshoots for the boundary streamlines. A comparison of the exact 
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solution with the computed results is shown in table 1. The incident shock 
angle is computed quite accurately; the reflected shock angle differs by 
8-10%. This error in the Prandtl-Meyer boundary correction is caused by the 
violation of the isentropic condition as required by Abbett's scheme. 

If the boundary method of Book, Boris and Hain (ref. 12), based on simple 
reflection principles, is followed for the above problem, no converged solu- 
tion can be obtained. If the ramp angle was reduced to tan“^ (1/20) the re- 
sults shown in figure 19 were obtained. Large overshoots appear in the 
streamlines near the boundary. The streamlines indicated by NY = 1 and 
NY = 15 are 1/2 cell from the surface. The surface pressures can be obtained 
by extrapolation. This is not done here since the results would be even worse 
than shown . 

All of the two-dimensional results were obtained with the Shasta FCT 
algorithm. The third order FCT was not stable with the splitting procedure 
used in this report. The reason is not clear, but it seems that the coeffi- 
cient of the fourth-order diffusion term of the modified equation, that is, 

[ (Ax‘*/24At)c^ (« ^ - 1)] has the wrong sign (it is negative for €2 < 1) for 
linear stability. It is conjectured that the nonlinear stability dominates 
the linear stability for the one-dimensional case but not for two-dimensional 
splitting procedures. 


V. CONCLUDING REMARKS 


A Lagrangian method has been developed to solve the unsteady inviscid 
gasdynamic equations by the FCT method. Unfortunately the boundary conditions 
as given by the original FCT method were not satisfactory and these were 
worked out for the special problems considered in this report. Shock waves 
and other flow discontinuities were captured accurately and monotonically . 

The method works well, but only with the inclusion of the proper bound- 
ary conditions. The results are much better than those obtained by the mixed 
Euler-Lagrange methods of the Los Alamos group. The purely Lagrangian 
approach used here avoids the problem of mesh generation. However, we have 
avoided highly sheared flows so that the method did not require any mesh re- 
adjustment to maintain stability and accuracy, which can be a problem with 
Lagrangian methods. 
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APPENDIX A 


DERIVATION OF GOVERNING EQUATIONS OF GAS DYNAMICS 
FOR AN ARBITRARILY MOVING COORDINATE SYSTEM 


We start with the Euler Equations 


U + F + G =0 (Al) 

t X y 


with 





— ~ 


— — 


p" 


pu 


pv 


pu 


pu'^ + p 


puv 

u = 


; F = 


; and G = 



pv 


puv 


pv^ + P 


pE 


(pE + p)u 


(pE + p)v 


— — ' 


_ _ 


^ aw 


where 

p - (r - i)p(e 

We now want to transform to another arbitrarily moving coordinate system, 
call it (a, b, t) where 


a = a(x, y, t) 
b = b(x, y, t) . 

T = t 

and 

1 = 3(x, y) = j 
D 3(a, b) " 


(A2) 


is the Jacobian of the transformation (a, b, t) to (x, y, t). 

Viviand's transformation technique (ref. 14) for a conservative 
system of equations applies and results in 

A(jii) + + + 0 (A3) 

Since we have a = a(x, y, t) and b = b(x, y, t) 


dt ” 3x 3t 3y 3t 9t 
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where we require that a, b, t be independent variables. 

Setting 9x/3t = Ujjj and 9y/3t = v^, i.e., the coordinate mesh veloc- 
ity, we obtain 


and , similarly. 


3a 

at '^m ^y 


ORIGINAL PAGE IS 
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— = -u b - V b 


Substituting into equation (A3) and simplifying^ 

^ - SnU) + Jay(G - v^U)] 


3 

3b 


Jb (F - u U) + J by (G 
x m 


- V U) = 0 

m J 


(A4) 


Now we need to invert the partial derivatives a , a , b , and b , 

X y X y 


- 

^x J 


-X. 


a = — r~ : and 

y J 


and setting x = t results in 


^ -t 


l^JU (u - u 


) -f JU (v - V ) 

raj m 


-y, 

^x ” T 


h = 

y J 




+ ^ I - (JU) (u - u^) ^ + JU (v - v^) 

oD 1 m J m 


^ \ 


0 (A5) 


where 


JU = 





Jp 


0 

Jpu 


^bp 

Jpv 

; c = 


JpE 


p V ■ 



_ - 
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and 


0 


p(-y u + X v) 

Si 3 

System (A5) are the first four equations of system (1). The remaining 
six equations are the kinematic relations describing the mesh velocities. 

We have previously defined 


3x 

3t ~ % 


and 


3t 


V 

m 


Differentiating these two with respect to a and b obtains the last six 
equations of system (1). 
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APPENDIX B 


f* 


THREE DIMENSIONAL LAGRANGIAN GASDYNAMICS EQUATIONS 


The system of equations is 


U + F +G^+H + S= 0 

a b c. 


(Bl) 


with 


U = 


" Jp 


~ 0 “ 

Jpu 


P<Vc - ’'cS' 

Jpv 



Jpw 

-► 

P‘Ve - Vb’ 

JpE 

; F = 

T M 

pJu 

X 


0 

y 


0 

2 


0 


-> 

G = 


"" 0 


0 


~ 0 “ 

p(y^2^ - y_z^) 

c a a c 


P'S'a'b - 


0 

p(x 2 - X 2 ) 

a c c a 


- ’‘a'b> 


0 



P'^a^b - *b>'a> 


0 

pjv« 

; H = 

pJw^ 

; and S = 

0 

0 


0 


-u 

0 


0 


-V 

0 


0 


-w 

_ — « 


_ 


— 


The nine other kinematic relations obtained by differentiating the last 
three equations of system (Bl) with respect to a, b, and c have not been 
included. They are required only to keep the system strongly conservative. 
The Jacobian of the transformation (a,b,c,t) to (x,y,z,t) is 
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j = (x,y,z) 

9(a,b,c) 


= X 


<.y, 


b c 


c b 




a c 




The superscript M again refers to the velocity components in the moving 
coordinate system (a,b,c). Thus, 


where 


and 


u = au+av + aw 

X Z 2 

M 

V = bu + bv + bw 

X y z 

w = cu + cv + cw 

X y 2 

‘ j{Vc - Vb} 

“y ■ I { - ’‘b^c } 

“z ■ J { S>’c - W } 

b ='7'|y^ -yz i 
X Jt-^ca acj 

b =-v-(xz -xz I 
y J ( a c c a / 

••z ■ 1 - Vc ) 

J {^a^b ^b*a } 

c = -7 1 ~ X 1 

y J ( b a a b J 

c = T I X y, - X y > 
z J ( a b b a I 
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TABLE 1.- COMPARISON OF EXACT SOLUTION WITH THE PRANDTL-MEYER 
AND NORMAL MOMENTUM BOUNDARY CORRECTION 
FOR THE REGULAR SHOCK WAVE REFLECTION 




^2 

^3 

(de*g) 

®2 

(deg) 

Exact 

2.2 

1.87 

1.47 

37.2 

46.6 

Prandtl- 

Meyer 

2.2 

1.8 

1.42 

37.2 

50.3 

Normal 

momentum 

2.2 

1.81 

1.4 

37.8 

51.6 
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LAGRANGIAN - COMPUTATIONAL PLANE 


Figure 1.- Transformation of physical to computational plane. 
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(a) Step 1 - Diffused transport solution, 51 points. 

Figure 3.- Pressure profiles in the Lagrangian formulation 
solved by the FCT method. 
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(b) Step 2 - Antidiffused solution with no flux limiter. 
Figure 3.- Continued. 


40 


t = 0 



I 1 1 1 1 1 I I I I I ^ 

0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 


X 


(c) Step 3 - Flux-limited antidiffused solution. 
Figure 3.- Concluded. 
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Figure 4.- Rectangular network of the third-order FCT algorithm. 
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Figure 5.- The one-dimensional Riemann problem. 
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t = 0.5545 
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'igure 6.- Riemann problem, Eulerian coordinates, solved numerically with 
the MacCormack scheme. Pressure as a function of x and time t, 

500 points, CEL = 1.0, and reflecting boundaries. 
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t= 0 

t = 0.2888 



(a) Pressure as a function of x and t. Nonequal masses in the fluid 
cells, 500 points, CFL = 0.9, and reflecting boundaries. 

Figure 7.- Riemann problem, Lagrangian coordinates, solved 
with MacCormack's scheme. 
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t = 0.2706 
t = 0.5464 



Figure 8.- Riemann problem, Lagrangian coordinates, solved with MaCormack’ 
scheme. Pressure as a function of x and t. Equal masses in the fluid 
cells, 500 points, CFL = 1.0, and reflecting boundaries. 
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Figure 10.- Riemann problem, Lagrangian coordinates, solved with the 
generalized Lax-Wendroff scheme with artificial dissipation. 
Pressure as a function of x and t. Nonequal masses in the fluid 
cells, 500 points, and reflecting boundaries. 
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Figure 13.- Riemann problem, Lagrangian coordinates, solved by the optimum 
generalized Lax-Wendroff scheme. Pressure as a function of x and t. 
Nonequal masses in the fluid cells, 51 points, and reflecting boundaries. 
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X 


Figure 14.- Riemann problem, Lagrangian coordinates, solved by the SHASTA 
1/8 FCT scheme. Pressure as a function of x and t. Nonequal masses 
in the fluid cells, 51 points, and reflecting boundaries. 
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Figure 15.- Riemann problem, Lagrangian coordinates (Boris' version - 
March 1976), solved by the third order FCT scheme. Pressure as a 
function of x and t. Nonequal masses in the fluid cells, 50 points, 
and reflecting boundaries. 
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Figure 16.- Regular oblique shock-wave reflection pattern - 

exact solution. 
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PRESSURE PROFILES ALONG 
INDICATED STREAMLINES 



Figure 17.- Regular oblique shock-wave relection. Prandtl-Meyer 
boundary correction, 1/5 ramp, Mach no. = 2.2. 
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Figure 18.- Regular oblique shock-wave reflection. Normal momentum 
boundary correction, 1/5 ramp, Mach no. = 2,2. 


57 
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Figure 19.- Regular oblique shock-wave reflection. Boundary conditions 
computed by Boris' procedure (1975), 1/20 ramp, Mach no. = 2.2. 
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